##########################################################################################

library(data.table)
library(optparse)
library(dplyr)
library(ggplot2)
library(ggsci)
library(RColorBrewer)
library(ggpubr)
library(ggsci)
library(patchwork)
library(ggrastr)

##########################################################################################

option_list <- list(
    make_option(c("--input_file"), type = "character") ,
    make_option(c("--njmu_file"), type = "character") ,
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    input_file <- "~/20220915_gastric_multiple/dna_combinePublic/public_ref/combine/MutationInfo.combine.addMolecularSubType.rmMIX.tsv"
    njmu_file <- "~/20220915_gastric_multiple/dna_combinePublic/baseTable/STAD_Info.addBurden.MSI_MSS.addCNVType.rmMIX.tsv"
    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/finalPlot/revise/age"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

input_file <- opt$input_file
njmu_file <- opt$njmu_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################

dat_burden <- data.frame(fread(input_file))
dat_njmu <- data.frame(fread(njmu_file))

###########################################################################################
## 去除年龄未知的人
dat_burden <- subset( dat_burden , !Age %in% c("[Not Available]" , "unknown") )
dat_burden$Age <- as.numeric(dat_burden$Age)
## 整合分子亚型
dat_burden$Molecular.subtype <- ifelse( dat_burden$Molecular.subtype %in% c("EBV" , "unknown") , "Unknown" , dat_burden$Molecular.subtype)
dat_burden$Molecular.subtype <- ifelse( dat_burden$Molecular.subtype %in% c("MSI" , "POLE") , "MSI/POLE" , dat_burden$Molecular.subtype)
## 整合MSS
dat_burden_mss <- subset(dat_burden , Molecular.subtype != "MSI/POLE")
dat_burden_mss$Molecular.subtype <- "MSS"

dat_combine <- rbind( dat_burden , dat_burden_mss )
dat_combine <- subset( dat_combine , Molecular.subtype!="Unknown" )

###########################################################################################
## njmu的提取IM的
dat_njmu <- subset(dat_njmu , Class=="IM")
## 注释IM的分子分型
tmp <- dat_njmu %>%
group_by(Patient , Age , Class ) %>%
summarize( BurdenExon = median(BurdenExon) )
tmp <- merge(tmp , dat_combine[,c("Tumor" , "Molecular.subtype")] , by.x = "Patient" , by.y = "Tumor")
colnames(tmp)[1] <- "Tumor"

dat_use <- rbind(dat_combine[,c("Tumor" , "Age", "Class" , "BurdenExon" , "Molecular.subtype" )] , tmp)

###########################################################################################
## 按照不同来源，对不同分子亚型对突变负荷的影响进行你和
y_lab <- "Mutation rate per MB"
plotAge <- function(dat_use = dat_use , molN = molN , y_lab = y_lab){
    plot <- ggplot(data = subset(dat_use , Molecular.subtype == molN) , mapping = aes(x = Age ,y = BurdenExon , color = Class))+
    geom_point_rast()+
    scale_color_aaas() +
    ggtitle(molN) +
    xlab(NULL) +
    labs(y = y_lab , x = "Age") +
    stat_smooth(method="lm",se=T)+
    stat_cor() +
    theme(panel.background = element_rect(fill = NA, colour = "black", size = 1),
        legend.position ='right',
        panel.grid.major = element_line(colour=NA),
        legend.text = element_text(size = 10,color="black",face='bold'),
        axis.text.y = element_text(size = 10,color="black",face='bold'),
        axis.text.x = element_text(size = 10,color="black",face='bold'),
        axis.title.y = element_text(size = 13,color="black",face='bold'),
        axis.title.x = element_text(size = 13,color="black",face='bold'),
        axis.line = element_line(size = 0.5))
    return(plot)
}

plot_mss <- plotAge(dat_use = dat_use , molN = "MSS" , y_lab = y_lab)
plot_gs <- plotAge(dat_use = dat_use , molN = "GS" , y_lab = y_lab)
plot_cin <- plotAge(dat_use = dat_use , molN = "CIN" , y_lab = y_lab)
plot_msi <- plotAge(dat_use = dat_use , molN = "MSI/POLE" , y_lab = y_lab)

#plot <- plot_mss + plot_gs + plot_cin + plot_msi + plot_layout(nrow = 4, byrow = FALSE)
plot <- plot_mss + plot_msi + plot_layout(nrow = 2, byrow = FALSE)
out_name <- paste0(out_path , "/MutBurden.Age.cor.pdf")
ggsave(file=out_name,plot=plot,width=6,height=4)



